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Summary 

Following the previous report, the proposed investigation on a Matched Asymptotic 
Expansion (MAE) method was carried out. It was concluded that the method of MAE is not 
applicable to launch vehicle ascent trajectory optimization due to a lack of a suitable stretched 
variable. More work was done on the earlier regular perturbation approach using a piecewise 
analytic zeroth order solution to generate a more accurate approximation. In the meantime, a 
singular perturbation approach using manifold theory is also under current investigated 

Work on a general computational environment based on the use of MACSYMA and the 
weak Hamiltonian finite element method continued during this period. This methodology is capable 
of the solution of a large class of optimal control problems. This part of the work continued until 
the departure of Dr. Robert R. Bless, who was supported under the grant as a Graduate Research 
Assistant at Georgia Tech. The first version of his computer code is now complete. A NASA 
contractor report (CR), based on Dr. Bless’ Ph.D. Dissertation [1] is presently in press. It contains 
the details of the general code as well as sample input and output. These details are not repeated 
herein. Work has continued since his departure to more fully understand the accuracy and 
limitations of the method and to adapt Dr. Bless’ code to use Mathematica which is available on a 
wider variety of computers than MACSYMA. 
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1. Research Accomplishments 


1.1 Matched Asymptotic Expansion (MAE) Investigation 

The MAE approach was first investigated to handle the launch vehicle atmospheric flight 
phase where the earlier regular perturbation approach did not produce satisfactory results. The 
essence of the MAE approach is outlined below. 

We first evaluate the outer solution which corresponds to a propulsion dominant phase of 
flight. This part of the solution is just our zeroth order solution in the earlier regular perturbation 
approach [12]. Next, we formulate the inner solution which corresponds to an aerodynamic force 
dominant phase. This part of the solution has been developed in [9] where the analytic solutions 
involve elliptical integrals of the first and second kind. Finally, a composite solution is formed by 
joining the outer and inner parts with the extraction of the common part using the matching 
conditions (see [10] and [1 1] for details). 


First of all, we non-dimensional ize and rewrite the dynamics as: 
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The notation of the variables are self-explanatory. The hatted variables are non-dimensional and 
the subscript i stands for initial value of the variables. Here 8 is a small physical parameter 
whereas 5 is a bookkeeping perturbation parameter with a nominal value of one. We are actually 
using a combination of singular and regular perturbation expansions. Setting £ to zero, we retrieve 
the zeroth order outer dynamics (no aerodynamic forces). On the other hand, introducing the 

✓s /\ 

stretched variables t = t/e,h = h/e and setting 8 to zero we obtain the zeroth order inner dynamics 
(no thrust terms). The atmospheric pressure and the Mach number dependency of the 
aerodynamics data will be introduced in the first order correction which will subsequently involve 
solving a set of non-homogeneous linear O. D. E's. The advantage of this approach over our 
earlier perturbation approach lies in the fact that we are now able to introduce aerodynamic forces 
in our zeroth order formulation. 


Our first attempt was to evaluate the composite zeroth order solution using the existent 
results in [9,12] by solving a set of 21 nonlinear algebraic equations. However, we were not able 
to find a solution The problem is not due to numerical difficulties but lies in the flaws of our MAE 
arguments. From the optimal solution using BNDSCO we determined that magnitude of the 
aerodynamic forces is less than 15% of the thrust over the whole trajectory, which means there is 
never a flight phase where the aerodynamic forces dominate over the propulsive force. However, 
the magnitude of the aerodynamics forces is largely determined by the dynamic pressure profile. 
The aerodynamic forces increases as dynamic pressure initially builds up due to gain in velocity. 
As the launch vehicle rises in altitude, the drop in air density outweighs the gain in velocity and the 
dynamic pressure decreases. This phenomenon indicates that we need two different regions (2 
pairs of outer and inner solutions) to formulate our whole trajectory (see figure 1). We also need a 
new independent variable such that if it is set to the right and left hand side limits, the two 
respective outer solutions are obtained. Clearly, altitude is not the suitable candidate because we 
can only retrieve the right hand side of our solution as h -> <». 
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In a nutshell, we conclude that the traditional (using altitude as the stretched variable) MAE 
approach which has found success in aero-assisted orbital transfer application is not applicable in a 
straight forward manner to the ascent trajectory launch vehicle problem. Further research is needed 
to identify a more suitable independent variable. Rather than pursue this line of investigation, we 
decided to return to our earlier regular perturbation study, and to investigate a singular perturbation 
approach based on a slow manifold concept. 

1.2 A Modified Regular Perturbation Approach 

An idea to extend the earlier regular perturbation method into the atmospheric flight phase is 
through a finite element approach. Since we cannot find a complete analytic zeroth order solution 
that incorporates the aerodynamic effects, our approach is to improve accuracy with minimal 
increase in computational complexity. Using several pieces of simple solutions instead one 
complete and complicated solution, we are able to improve the zeroth order solution so that it 
accounts for the aerodynamic effects. 


From our earlier study, the state dynamics can be fairly represented by those of a flat Earth 
no atmosphere approximation. However, this is not true of the costate dynamics. If we use the 
previous approximation, we will end up with (in a rectangular coordinate system): 
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The approximate set of zeroth order dynamics are especially poor in X u and Xj because both 
derivatives become zero to zeroth order in e. Consequently, they produce large forcing function 
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terms (in L 2 -norm sense) in the first order correction dynamics which may cause divergence of our 
corrected results. The easy way to decrease these large error magnitudes is to represent the Xj, and 


Xj with linear function such that 
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and p u , p r are constants to be determined by other means. The optimal control of the zeroth order 
problem is now governed by a bilinear tangent law. 


The constants p u and p r are evaluated using collocation method [13]. We approximate the 
solution with pre-specified functions, in this case first order polynomials. Constraining the 
solution such that it is continuous at the node and satisfies the original dynamics at the mid point of 
each segment determines the unknown coefficients of the polynomials. Mathematically, these 
constraints are formulated as follows: 
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The unknowns to be solved are the nodal values (subscript 1, 2) of the interpolated variables. 
However, for this linear function case, we simply equate the unknowns p u and p r with the 
averages, ie. 
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Figures 2-8 show an open loop 4-piece zeroth order solution segmented at 30s, 60s, 90s. In our 
present formulation, we also treat X v as a linear function with the unknown p v . The first 3 
segments are computed using collocation method described above, and the last segment uses a flat 
Earth approximation (p u = p r = 0). Asa comparison, the costate profiles (figures 6 - 8) of the 
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earlier 1 -piece zeroth order solution and the optimal solution (see Errata) are also plotted. We can 
clearly observe the significant improvements of this modification. 


The new analytic expressions of v, u, h are given below: 
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The state transition matrix can be found by differentiating the above analytic expressions with 
respect to the initial values (c v , c u , c T are the initial costate values). First order correction using the 
zeroth order solution above and the state transition matrix will be obtained in next progress report. 


1.3 Singular Perturbation Approach Using Manifold Theory 

In [15], we showed that a singular perturbation, using a 2-state model with mass and 
energy as slow variables, failed because the flight path angle dynamics are highly coupled with the 
slow variables at high flight path angle. However, if we use a more accurate 3-state model (mass. 
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energy and altitude), a chattering solution of flight path angle will be encountered. Our proposed 
research for the next reporting period is to attempt to use the Manifold condition [14]: 

g = £ (|*f+|*> 20) 

ox at 

where 

— = f(x,z,t) e^ = g(x,z,t) z = <{>(£, x,t) 

dt dt 

to generate another zeroth order outer solution (slow manifold). Since we now include e in our 
slow manifold (<|>) formulation, the chattering effect is eliminated. Our first step is to demonstrate 
that we can compute an off-line slow manifold solution and perform an on-line boundary layer 
(inner solution) correction for the flight path angle dynamics. This will result in a nonlinear 
feedback control solution for the angle of attack (see below). 

H = XEoE^o.ho.mo.aJ + XhoVoSiny-Xmok + X^Eo.ho.mo.a) = 0 (21) 

H a =0 (22) 

These two equations are used to determine a and The subscript 'o’ stands for the initial value. 

1.4 Finite Element Analysis 

The main accomplishment during this reporting period involved the development of the 
general code. The main purpose of the general code is to reliably solve a large class of optimal 
control problems with a minimum of user-written subroutines. To this end, the general code runs 
on a SUN 3 and later workstations. It and requires a FORTRAN 77 compiler, MACSYMA [2], 
and the Harwell subroutine library [3]. The general procedure can be broken into three parts that 
must interface together. The first part is the FORTRAN code. This code contains all the 
subroutines necessary to solve any of the optimal control problems described above. However, if 
certain problems require table look-up routines (such as aerodynamic data for a rocket model), then 
these subroutines must be given by the user and interfaced to the rest of the general code. Thus, 
there may be a need for some user programming for certain problems. The second part of the 
general procedure is the use of MACSYMA. The user must supply an input file specifying the 
problem. This input file is in symbolic form and will be loaded into MACSYMA. MACSYMA 
will then evaluate all the necessary expressions and automatically generate the FORTRAN code. 
This code is spliced into a template file and becomes one of the subroutines. The third and final 
part of the general procedure will consists of subroutines to generate initial guesses that will 
reliably converge. A continuation method is being used which converts the algebraic equations to 


8 


initial-value ordinary differential equations. A second-order Runge-Kutta method is then used to 
integrate the equations and obtain initial guesses for a Newton-Raphson method. We also 
continued to further document the methodology through the publication of one paper [4] and the 
completion of three others. The first of these three is a technical note [5] which covers the 
extension of the method to state-control inequality constraints. The second deals with the 
application of the method to the ALS problem, per se [6]. Both of these are now accepted for 
publication. The third deals with the general code and will be presented at the upcoming ACC 
meeting [7]. 
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2. Errata 


An error in our previous 2-stage ALS optimal solution was found. The optimality 
condition was incorrectly formulated due to a missing conversion factor from degree to radian. 
The correct results are now shown in Figures 9-15. There are 3 jumps in the control profile 
(Figure 9). The first two jumps are due to the fact that the Hamiltonian is a non-convex function of 
the control. These jumps occur at about Mach 1 .4 and 2.0 respectively. The last small jump is due 
to staging. However, the costates are all continuous. Though the control profile changes 
substantially, the performance index (final time) differs by less than 0.1s. Figures 16, 17 are the 
optimal solution under aq-constraint. The final time in this case is 362.103s which is 0.007s more 
than the unconstrained case. 
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3. Future Research 


3.1 Analytical Investigation 

We will follow two different directions. One will be the continuation of the modification of 
the regular perturbation technique. A first order closed loop simulation will be done. We will 
investigate the effect of number of segments and the segment intervals on the computational 
performance. At present a first order correction can be done in 0.5 to 7.0 CPUsec, depending on 
the current vehicle altitude on a MacIIci (a 32-bit 25MHz PC). At low altitude, the complicated 
aerodynamic effects require a more dense integration steps to complete the quadratures, however, 
the computational time can be significantly shortened if we perform parallel processing. The other 
line of research direction is to investigate the slow manifold approach to singular perturbation to a 
launch vehicle trajectory. Further approximation and analysis are expected. 

3.2 Finite Elements Work 

In the finite element analysis area we plan to continue to port the general code and complete 
a user's manual for it. We further plan to document the methodology through the completion of 
one paper (which we are now revising in response to reviewers) on the application of the method 
to launch vehicle trajectory analysis, two technical notes on control and state inequality constraints, 
one paper on the general code, and a user’s manual for the code. We continue to receive calls from 
parties interested in application of the methodology in industry, and still hope to transfer the 
technology to an industry application in the future. 

In a future phase of the work we hope to extend the work to higher-order finite element 
shape functions - the so-called p-version of the finite element method. This approach has been 
shown to be of value in linear time-domain problems [8] but has not yet been investigated for the 
nonlinear case. 
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L. H. S. Outer Solution 
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